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ABSTRACT 

We study the role resonant scattering plays in the transport of Lya photons 
in accreting protoplanetary disk systems subject to varying degrees of dust set- 
tling. While the intrinsic stellar FUV spectrum of accreting T Tauri systems 
may already be dominated by a strong, broad Lya line (~80% of the FUV lu- 
minosity), we find that resonant scattering further enhances the Lya density in 
the deep molecular layers of the disk. Lya is scattered downwards efficiently 
by the photodissociated atomic hydrogen layer that exists above the molecular 
disk. In contrast, FUV-continuum photons pass unimpeded through the pho- 
todissociation layer, and (forward-)scatter inefficiently off dust grains. Using 
detailed, adaptive grid Monte Carlo radiative transfer simulations we show that 
the resulting Lya/FUV-continuum photon density ratio is strongly stratified; 
FUV-continuum dominated in the photodissociation layer, and Lya-dominated 
field in the molecular disk. The enhancement is greatest in the interior of the 
disk (r ~ lAU) but is also observed in the outer disk (r ~ lOOAU). The majority 
of the total disk mass is shown to be increasingly Lya dominated as dust settles 
towards the midplane. 

Subject headings: Stars: variables: T Tauri, Herbig Ae/Be - Protoplanetary disks - 
Methods: numerical - Radiative transfer 
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Introduction 



T Tauri stars (TTS) are pre-main sequence low mass stars frequently described as 
young analogs of our Solar System. Surrounded by circumstellar disks of gas and dust, 
these systems mark the earliest stages of p l anet-formation, a process now understood to b e 



common in our Galaxy (jAdams et al. 



1987 



Kenyon fc Hartmann 



1995 



Pollack et al. 



19961). 



Due to their close proximity to the parent star, intense ultraviolet and X-ray radiation 
fields greatly impact the evolution of the protoplanetary disk. The far-ultraviolet (FUV, 
hu < 13.6eV) field is of particularly bro ad interest as it provides the rmodynamical gas 



heating through the photoelectric effect (jWeingartner &: Draine 



200 ll ). and affects disk 



chemistry directly through the photo-dissociation and -ionization of molecules and atoms 



( JBergin et al. 



20071 . and references therein). FUV photodesorption of ices contributes to 



the transport of material bet ween solid and gas phases, whilst enabling complex molecule 



formation on grains surfaces (jOberg et al. 



2009al Jbl). Alongside X-rays and cosmic-rays, the 



FUV field is also involved in the ionization balance of the disk; of central importance in 
dynamical proces ses such as the magneto-rotational inst ability that is believed to play a 



role in accretion (JBalbus fc Hawley 



1991 : 



Gam 



mie 



19961). Eventually the dispersal of disk 



20061 ). although resolving 



gas may be driven by UV photo-evaporation ( 1 Alexander et al. 

the contri butions made by the separate FUV, EUV an d X-ray components is still a matter 



of debate (IGorti fc HoUenbach 



2009 



Owen et al.lboioh 



One feature of the FUV spectra in T Tauri systems which separates the accretors from 
the non-accretors, is a strong, broad Lya line (accounting for as much as 80% of the total 
FUV lur ninosity) . In riearby systems such as TW Hydra the Lya line can be detected 



directly (JHerczeg et al. 



20021 ). whilst in obscured systems it s 



through the excitation of H2 molecules ( jHerczeg et al. 



200 



presence may be inferred 



Bergin et al. 



20041 ). It was 



soon realized that the concentration of UV photons in a relatively narrow wavelength range 
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'1215.67 ± 1 A) would have important implications for wavelength- dep ende nt chemical 



20031 ). For example, 



processes, specifically the photodissociation of molecules (iBergin et al. 

the molecules HCN, OH, H2O and H2CO can all be photodissociated by a Lya photon, 

unli ke the closely related spe cies CN, 0H+, H2O"'" and HCO"'' (for a more complete list. 



see 



van Dishoeck et al. 



20061 ). Abundance ratios frequently used as chemical diagnostics. 



such as A^('CN)/iV (HCN), will be cha nged (in this case increased) by Lya-dominated 



photodissociation. 



Fogel et al 



( 1201 ll ) have shown quantitatively that the inclusion of 



differential photodissociation by Lya ca n some instance s help bring abundance ratios 



into better agreement with observation (JThi et al 



20041 ). While the the presence of Lya 



decreases the abundance of HCN, NH3, and CH4 by an order of magnitude or more, other 
species such as CO2 and SO are enhanced. SO actually has a significant cross-section 
at 1216A, however so does SO2, and it is the photodissociation of SO2 which partially 
replenishes the SO abundance. In this instance, the presence of Lya actually increases the 



Fogel et al 



(1201 ih models 



SO abundance. Similarly, the abundances of H2O and OH in the 
increased with the introduction of Lya, even though both H2O and OH have a significant 
photodissociation cross-sections at 1216A. In this case the gas phase abundances are more 
than adequately replenished by efficient Lya photodesorption from grain surfaces. Further 
chemical modeling is required before the implications of Lya-dominated photochemistry 
are fully understood. 

Despite the interesting chemical implications of a Lya dominated FUV field, the 
problem of precisely how it penetrates the disk has received little attention in the literature. 
In this paper we attempt to shed light upon its basic transport. Although frequently fiared, 
protoplanetary disks remain relatively fiat objects with vertical s calehe ights that are small 



compared to the radial extent of the disk 



Kenyon &: HartmannI (119871 ). As a result, the 



majority of the mass (residing near to the disk midplane) is concealed from the parent 
star. Visual optical depths at the midplane, measured along lines directly from the star. 
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typically exceed 10^, implying that essentially no stellar photons can penetrate these parts 
directly. Instead, radiatio n must be scattered downwards from the surface of the disk 



d 



van Zadelhoff et al 



20031 ) ■ Scattering is therefore of central importance when discussing 



the transport of stellar radiation in protoplanetary disks. 

The penetration and absorption of X-ray and FU V radiation determ ines much of the 



Aresu et al 



2011 



and references 



thermodynamical stratification of the upper disk (e.g. 
therein). For simplicity we can identify three layers; the uppermost layer being a 
photodissociation layer illuminated directly by the intense, unattenuated stellar radiation 
field - this layer is dominated by atoms, since the photodissociation timescale of molecules 
is very short. As unscattered photons penetrate into the disk they eventually reach an 
irradiation surface, defined to be where the optical depth integrated from the star is of 



order unity, r* ^ 1 (ICalvet et al. 



1991 



Chiang &: Goldreich 



1997 



Watanabe fc Lin 



20081 ). 



It is here that, in an average sense, photons undergo their first interactions with the 
disk. In the case of interactions with dust grains, this results in attenuation through both 
absorption and scattering; the latter leading to the emergence of a diffuse component to 
the radiation field. Since the radiation field decreases with depth, while the local mass 
density increases, we eventually enter an intermediate warm molecular layer, where H2, 
CO and a host of other gas-phase molecules may exist in a dynamic equilibrium with the 
photodissociating FUV field. It is in this warm molecular layer that the photochemical 
processes are most important. Despite its attenuation, the radiation field is sufficient to 
maintain grain temperatures above 30K, thus preventing widespread freeze-out of molecules 
onto grain surfaces. At even greater depths we enter the optically thick, dense midplane 
region containing > 99% of the disk mass. Here the stellar FUV field is highly attenuated 
and consists entirely of scattered photons. With the exception of the innermost few AU, 
the tem peratures may be low enough (< 30K) to permit the freeze-out of several molecular 



species (JAikawa et al. 



2002 



Bergin et al 



20071 ) 
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The transport of FUV photons is further comphcated by the gradual 



settling of dust grains as the disk evolves towards planetesimal formation (JThroop et al. 



growth and 



2001 



Wilner et al. 



20051 ) ■ Not only does grain-growth affect the optical properties of 



individual grains, but the radial and ver tical drift of so 



segre gation of the grain-size population (IWeidenschilling 



ids relative to gas causes a 



1977 



DuUemond fc Dominik 



20041 ) ■ The majority of the grain opacity at FUV wavelengths is attributable to relatively 
small grains with radii a < l/xm, which subsequently reemit absorbed energy at infrared 



wavelengths, making it possible to estimate their abundance (JDominik fc DuUemond 



20081). 



Observationally, the abundance of small grains (per H nucleus) relative to that found in the 
ISM, e, typically takes a value much less than one; for example, in the Taurus star-forming 



region the median value is e ~ 0.01 (JFurlan et al. 



20061). 



Our basic picture of FUV transport is somewhat incomplete once we include Lya 
photons. In view of the large cross-section due to resonant scattering by H atoms, cr-^Yai 
it seems entirely plausible that Lya photons experience their first interactions not with 
dust grains but with H atoms residing in the uppermost photodissociation layer of the disk 
(hereafter 'H-layer'). The irradiation surface for Lya photons is therefore defined as the 
location where r* = j crj^^^n^H) ^ 1. The precise location depends quite sensitively on the 
width of the stellar Lya line profile, although the Lya photons that propagate freely to the 
disk surface are generally those in the wings of the scattering line profile presented by the 



disk ( iHerczeg et al. 



20021 ) . Even hundreds of Doppler widths out in the Lorentz wings of the 



scattering line profile the cross-section may still be significantly larger than that due to the 
(possibly settled) dust population. Figure [T] shows the Lya line from TW Hya, plotted next 
to an estimate of the scattering line profile due to warm H atoms in the disk. Assuming the 
photodissociation layer supports H columns of N*{H) > 10^°cm~^ (integrated along a line 
from the star), it seems probable that the Lya field will first be scattered by this atomic 
H-layer. The optical depth of dust in this layer will be smaller by a factor of approximately 
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0.1e« 1. 

In stark contrast to Lya transport, FUV-continuum photons will stream unimpeded 
through the H-layer, eventually striking the disk at the conventional irradiation surface 
defined by dust. Here they will begin to scatter, however scatteri ng by typical i nterstellar 



2004h . Multiple 



dust populations is considered to be somewhat forward-throwing (iGordonl 

scatterings are therefore required before an appreciably diffuse component is generated, and 

accompanying this scattering will be a degree of pure absorption. 

In this paper we explore numerically the basic phenomenology of resonantly scattered 
Lya in disk models. An emphasis is placed on contrasting the transport of Lya with that of 
the FUV-continuum field in the vicinity of the Lya line. We proceed with Section [2] where 
we discuss the properties of our disk models. Radiation transfer of Lya and FUV-continuum 
photons is dealt with in Section [31 This includes an approximate - but largely self-consistent 
-computation of the H distribution necessary for Lya transport. Results are presented in 
Section HI The paper concludes with a discussion and summary. 

2. Disk models 



We base our disk structures on a suite of three disk models from 



D'Alessio et al. 



(120061). The disks are differentiated by the depletion of small grains in their upper layers, 
e = {0.01,0.1, 1.0}. Physically, the parameter e represents the ratio of dust concentration 
relative to that found in the interstellar medium. Although the models include the effects 
of dust depletion on the hydrostatic disk structure, they lack a separate thermodynamic 
treatment of the gas, which becomes thermally decoupled at the low densities found 
in the disk atmosphere. While the D'Alessio et al. disks suffice for illustrating the 
essential features of Lya and FUV-continuum transport, a more self-consistent treatment 



will eventually require the inclusion of detailed coupling between the radiative transfer, 
thermodynamics and disk structure. Several disk models exist that already include the 



thermal decoupling o f gas and dust that 



Ercolano et al. 



2009 : 



Gorti fc HoUenbach 



generally leads to 



2009 



Woitke et al. 



arger 



2009; 



disk scaleheights (e.g. 



Owen et al. 



2010|). We 



assign a nominal FUV luminosity of O.OIL© to the central star and a kinetic temperature 
of Tg = lOOOK to the gas in the upper layers of the disk. In view of these modifications 
we refer to the models as 'modified D'Alessio disks'. A cross-section of the total hydrogen 
density, nn = n(H) + n{B.2), in one of the disk models is shown in Figure [2j 

The grain-size distribution used in the original D'Alessio calculations follows 



Mathis et al. 



(119771 ). dn/da oc a"^-^, where a is the grain radius. In fact there are two 
populations of grains ('large' and 'small') in the D'Alessio models, however the large grains 
are confined to the disk midplane and play no role in what follows. For the 'small' grain 



population the grain-size limits are a^ 



O.Olyum and amn-r = 0.25Mni 



These limits are consistent with those proposed by 



grams. 



Mathis et al. 



'double check this). 



(Il977l ) for interstellar 



Relevant optical properties of the 'small' grain population are shown in Figure [31 In 
this paper we are primarily concerned with FUV-continuum wavelengths in the vicinity 
of 1215. 67A. At this wavelength the 'small' grain population exhibits an appreciably 
forward-throwing phase function (a symmetry factor of q ~ 0-6) • The wi dely used 



Henyey-Greenstein phase functions (JHenyey fc Greenstein 



1941 



Witt 



19771 ) evaluated for 



g = 0, 0.4 and 0.6 are shown in Figure Si 

The computational domain used for the subsequent radiative transfer simulations 
comprises the modified D'Alessio disks immersed in a background mesh. The modified 
D'Alessio disks are spatially 1+lD datasets; the purpose of the background grid of nodes is 
to transform these models into true 2D distributions. The resulting discretization consists 
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of a list of nodes, each of which has associated with it a location as well as relevant physical 
quantities (e.g. densities) which are interpolated from the original D'Alessio models. The 
background mesh supports unstructured distributions of nodes in order to accommodate 
spatially irregular datasets such as the D'Alessio models. To the extent possible, the initial 
background grid exhibits grid refinement such that its spatial resolution matches that of 
the D'Alessio data. 

The nodes (constrained to lie in the r — z plane in a cylindrical coordinate system) 
ultimately form the vertices of triangular cells. Figures E] and E] give the reader an impression 
of the spatial distribution of nodes and how they result in a tessellation of triangular cells. 
The triangular cells are a result of the connectivity obtained by performing a Delaunay 



tessellation on the node list 



can tessellate the volume (jOkabe et al. 



. In 2D the tri angle is the simplex: the simplest shape that 



19921 ). In 3D the simplex is the tetrahedron. There 
are many alternative ways to connect the nodes that give rise to a sensible tessellation of 
triangular cells, however the Delaunay tessellation has unique properties that arguably yield 
the optimal connectivity. First, the connections link those nodes that compete for space. 
This is most easily understood by considering the corresponding Voronoi diagram (the 
dual of the Delaunay tessellation). Second, the resulting triangles maximize the minimum 
internal angle, thus reducing the frequency of skinny triangles. While the meshes in Figs. 
|5] and |6] ostensibly appear regular, they are treated as fundamentally unstructured. Had 
the background mesh been cast randomly (e.g. according to a Poisson point process) it 
would be possi ble to shed all traces o f the underlying coordinate system in the spatial 



discretization fiRitzerveld fc Icke 



20061 ). For these reasons the Delaunay tessellation is 
becoming an attractive choice for generating high-quality, unstructured discretizations of 
spatial objects. 



^See, for example, http://www.qhull.org/ 
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3. Radiative transfer 

Very few radiative transfer problems of astr ophysical interest can be solved analytically 



1960l ). The transport of 



without resorting to significant simplifications (I Chandr asekharl 
radiation in protoplanetary disks is no exception. In view of the geometrical complexity 
and assortment of opacity sources, we opt to solve the transport problem using the highly 
versatile Monte Carlo simulation method. Numerical details of our Monte Carlo radiative 
transfer code are described in the Appendices. In this section we restrict the discussion to 
aspects which are relevant to the physical context. 

3.1. Transport of FUV-continuum 

The transport of FUV-continuum photons is assumed to be governed solely by 
interactions with dust grains. At extreme levels of dust settling (e < 0.001) it is possible 
that Rayleigh scattering by II2 molecules becomes relatively important. Anisotropic 
scattering by grains is treated using the one-parameter Henyey-Greenstein phase function, 
the angular scattering probability distribution of which is given by. 



m ' '-'' 



A7T[l + g^-2gcosef/^' ' ' 

where g = [—1, 1] is the anisotropy parameter. Isotropic scattering corresponds to g = 0, 
forwards scattering to g = 1, and backwards scattering to g = —1 (Figure H]). The 
single-scattering albedo cu ~ 0.4 indicates that these grains have a slight preference towards 
absorption rather than scattering. From u and g we can conclude that multiple scatterings 
will be required in order to generate a diffuse field that can penetrate the disk vertically. 
Accompanying these scatterings will be significant absorption. 
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3.2. Transport of Lya. The H/H2 distribution 

The transport of Lya is treated like FUV-continuum but with the addition of resonant 
scattering. Accordingly, it is necessary to first determine the distribution of atomic 



hydrogen. The formation of 



of dust grains (ISavage et al. 



iy, in dense, dusty environments usua 



1977 



Spitzer 



1978 



Cuppen fc Herbst 



y occ urs on the surfaces 



20051 ). The destruction 



of H2 occurs with ~ 10% probability following an absorption of a FUV p hoton into 



HoUenbach fc Tielens 



I999I). Since the 



the Lyman- Werner band system (912 — llOOA, 

photodissociation of H2 occurs through line absorptions, these transitions become optically 

thick over very modest columns, A^(H2) ~ lO^^cm"^. Beyond this point we see a run-away 



formation of 'self-shielded' H2. While any photo-destruction do minated process that 



competes with formation can in principle result in self- shielding (iBethell fc Berg: 



20091) 



the concentration of H2 opacity in discrete lines makes it readily self-shi elding. Balancing 



formation and destruction rates gives the steady-state molecular fraction (ISpaans fc Neufeld 
1997al ). 



f{H,) ^ 



n(H2 



n-^Rn^e 



n 



H C + 2ngi?H2e 



(2) 



wher e Rh2 is the high temperature H2 formation rate coefficient (jCazaux fc Tielens 



2004 



. For the self-shie l ding 



20101), and ( is the H2 photodissociation rate ( = (qFs[N (H.2) 
function, Fs[A^(H2)], we adopt the closed-form expression from 
although c ompa rison simulations were also performed using the tabulated results of 



Draine fc Bertoldil fll996h 



Lee et al 



(119961 ). These approximations were shown to result in H-layers with comparable 
geometries and column densities. In the context of our Monte Carlo simulation we 
eval uate the self-shielding fun ction along the length of each photon packet trajectory 



[e.g. 



Spaans fc Neufeld 



1997bl ). The inclusion of a full treatment of H2 photodissociation 



i s currently 
(JRolhg et al 



peyon d the capabilities of most multidimensional photochemical codes 
2003). 
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3.2.1. Iterative H/H2 calculation 

The calculation of the H/H2 distribution is an iterative one, since to compute n{Y{2) 
at a point requires us to know the degree of self-shielding due to n(H2) everywhere 
else. An initial distribution of H and H2 is required to start the iterative calculation. 
From this our H2 will 'grow' as self-shielding is established over the course of repeated 
iterations. Ideally, this initial condition should be constructed to underestimate the H2 
density: if H2 is overestimated then regions that are artificially self-shielded will be slowly 
photodissociated away over the course of the iterative process. In such over-shielded cases, 
a layer of excess H2 with a thickness corresponding approximately to the self-shielding 
column (A^(H2) ~ lO^'^cm"^) will be photodissociated in each iteration. If the excess H2 
is distributed over a large volume, the convergence towards the true H2 distribution will 
be very slow. Since the introduction of self-shielding always increases the H2 density, we 
can set a strict lower limit by computing n(H2) using Eqn |2] ignoring all forms of shielding 
(this includes dust attenuation, although we retain the geometrical inverse- square dilution 
of radiation). The resulting n(H) is shown in the left-hand panel in Figure |5l Even without 
self-sheilding, atomic hydrogen is limited to the upper disk layers due solely to the increase 
in density that promotes H2 formation deep in the disk. 

The iterative process commences with a full radiative transfer simulation at a single 
representative Lyman- Werner band wavelength (centered at A = lOOOA). The accumulation 
of II2 column density along photon packet trajectories is recorded, and information 
regarding the self-shielded radiation field is deposited in each cell traversed by the photon 
packet. Dust scattering and absorption are both included. However, since H2 shielding 
proceeds so vigorously, dust shielding plays a relatively minor role in determining the final 
II/H2 distribution. This is especially true in dust-settled scenarios, although a rigorous 
consideration must establish the criteria that separates purely H2 self- shielding situations 
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from those in which H2 self-shielding is initiated by dust attenuation ( IWolfire et al.lll995[ ). 



Once approximately 10^ photon packets have been run, the resulting self-shielded 
radiation field is used to compute a new n(B.2) using Eqn |2l Since the original grid is 
generally too coarse to properly resolve the emerging H/H2 transition we refine the grid 
before proceeding with a new iteration. If the H2 distribution is described on grid that is 
too coarse, individual cells may self-shield, causing a local overestimation of n{B.2) that 
may subsequently propagate through the grid. The refinement criterion is based upon the 
fractional change in n{B.2) between iterations. This ensures that the spatial resolution of 
the grid is increased in regions which are most sensitive to changes in n(H2). Regions that 
are either entirely unshielded or totally self-shielded (?2(H2) — )■ 0.5?2j|) will not receive grid 
refinement. The refinement of a cell is achieved by placing a single node in the geometric 
center of the cell. Once new nodes have been added a Delaunay tessellation is performed on 
the augmented node list, resulting a new set of cells. Physical quantities are interpolated 
onto this new grid and the next iteration cycle is begun. 

The criteria for a converged solution must not only require a solution that does not 
change significantly with further iterations, but also that every cell in the H/H2 transition 
must not contain a sufficiently large A^(H2) to self-shield all by itself. For a cell to be 
considered part of the H/H2 transition it must exhibit a n{B.2) that is greater than its 
unshielded value and less than its maximal possible value of O.Snjj. 

The nodes and n{}i) distributions after the first and second iterations in the e = 0.01 
model are shown in Figure O It is worth noting that the uppermost layers are entirely 
photodissociated so that n(II) ^ n-^ throughout the iterative procedure. It is the lower 
boundary of the H layer that is affected by self-shielding, and this can be seen to rise in 
height as the bulk of the disk gradually becomes molecular. For the modified D'Alessio 
disks we typically undertake 15-20 iterations. Convergence is assumed when n{B.2) changes 
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by less than 5% everywhere in the domain. It is important to reahze that Monte Carlo 
noise will generate variance between iterations, even after many iterations. If the Monte 
Carlo noise level is known then Eqn|2]can be used to estimate its effect on n(H2). 

To summarize, the basic steps of the iterative H2 calculation are as follows; 

1. Compute an initial distribution (lower limit) of H2 due to dissociation by the 
unattenuated stellar field (i.e. ignoring H2 self-shielding and dust attenuation). 

2. Propagate photon packets using Monte Carlo methods, following H2 self- shielding and 
dust-absorption along trajectories. 

3. Compute new H2 distribution due to this photodissociation field. 

4. Add nodes to regions where the II/H2 transition is detected. Interpolate physical 
properties onto the new mesh. 

5. Repeat from step 2 until solution has converged. 



3.2.2. Lya transport 

Once the H-layer has been calculated it is possible to propagate Lya photons. The 
resonant scattering process is simplified somewhat by ass ociati ng with it a Voigt line profi 



x) , and an isotropic phase function. 



Ahn et al. 



fcoOlh and 



Laursen fc Sommer-Larsen 



(I2OO7I ) discuss the treatment of resonant n = 1 — )■ 2 transitions in more detail, although 
for our purposes the simplified treatment is sufficient q As discussed previously, most Lya 
photons from the TW Hya spectrum scatter in the Lorentzian wings of the Voigt profile (see 
Figure [1]). Unlike the Doppler core, which is relatively temperature sensitive, the wings of 



^A more complete microphysical treatment is required if polarization is to included. 
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the scattering profile do not change greatly with the temperature of the gas. This somewhat 
mitigates the ad hoc temperature imposed upon the gas in our modified D'Alessio disks. 

The convolution of the Gaussian (thermal broadening) and Lorentz (natural) line 
profiles gives rise to the Voigt profile, 



a{a, x) = f: 



e\/7r 



12" 



-H{a,x), 



(3) 



' mecAun 
where fu is the oscillator strength for the the 1 — ?■ 2 transition, e is the electron charge, c 

the speed of light, and nie the electron mass. Aud = v^vth/c is the thermal Doppler width 

where Vth is the mean thermal velocity of the scatterers and vq is the frequency at line 

center. H{a,x) is the Voigt profile defined by the integral. 



H{a,x) = — 

71 



+ 00 



exp —y 



-dy. 



(4) 



, (x + yy + a^ 

While the scattering is coherent in the frame of the scattering atom, an overall 
frequency change occurs due to the difference between the incoming and outgoing Doppler 
shifts required to change betw een atom and disk frames. This is known as partial frequency 



redistribution (Hummer 



19621 1. The fractional change in photon frequency is of order 
the thermal Doppler width Az/£) = v^Vth/c where Vth is the mean thermal velocity of the 
scatterers and z/q is the frequency at line center. Partial redistribution of the scattered 
Lya photon is include d in our code using microscopi c Mon te Carlo techniques based on 



Avery fc House! (119681 ) and IZheng fc Miralda-Escude (120021 ). However, since the intrinsic 



stellar Lya line shape is already very broad (~ 300 Doppler widths, see Fig. [T]), and the 
H-layer not extremely optically thick, little appreciable frequency redistribution of the Lya 
line is expected. As such one can largely ignore frequency redistribution in our models, and 
instead treat the resonant scattering as elastic. This simplification greatly expedites the 
calculation. We ignore systematic Doppler shifts due to the azimuthal velocity component 
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of the disk (^ 30kms"^ at lAU, 



Gayley et al. 



200l[ ) since for much of the disk this induces 



Doppler shifts that are modest compared to the stellar Lya line width (> SOOkms ^] 



4. Results 

The vertical distribution of hydrogen in the three disk models is shown in Figure [HI 
The molecular hydrogen column density A^(H2) shows the characteristic S-shape associated 
with the onset of self-shielding as we descend into the disk. The vertical column density of 
atomic hydrogen clearly varies with the degree of dust depletion, A^(H) ~ 10^^ — lO^^cm"^ 
- this is the column thickness of the H-layer that scatters Lya photons. Thicker H-layers 
are associated with dust depleted models, due to the proportional reduction in the H2 
formation rate. The precipitous drop in n(H) seen in Figure [8] is indicative of the sudden 
onset of H2 self-shielding. The vertical optical depth (for resonant scattering) of the H-layer 
depends upon the wavelength displacement from line center; for the Lya spectrum shown 
in Figure [1] the photon-averaged optical depths are Tt ~ 0.2,1,6 for e = 1,0.1,0.01, 
respectively. Although these vertical optical depths are not large, when viewed from the 
star the slant optical depth this H-layer is typically Tt ~ 20ry ^ 1 and thus the 
H-layer will intercept essentially every stellar Lya photon incident upon it. We are mostly 
interested in following the portion of Lya emerging from the lower surface of the H-layer. 
The remainder is scattered into space from the upper surface. It is worth noting that the 
H-layer is not so optically thick to these line- win g photons that they become line-trapped 



and destroyed by dust absorption (JNeufeld 



1990) • While these results suggest that resonant 



scattering of Lya is indeed important, the H-layers computed lack the self-consistency 
required to render them entirely realistic. 

Figure [9] shows spatial maps providing a visual comparison of the Lya and FUV- 
continuum photon densities in the e = 0.01 modified D'Alessio disk. In addition to the 
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photon density Uph we also show flux arrows (direction only), and the anisotropy of the field 
7 (the ratio of net flux to photon density). The highly flared H-layer and its isotropizing 
effect on the Lya field is clearly seen in the left-hand anisotropy panel. The net flux of 
transmitted Lya photons emerges from the H-layer perpendicularly, providing a more direct 
illumination of the disk below. 

Vertical profiles of the Lya/FUV-continuum photon density ratio are shown in Figure 
[To] for e = 0.01, 0.1, and 1, at radii r = 1 and lOOAU. The same qualitative behavior of the 
ratio is seen at all radii, although the quantitative effect is greatest in the inner disk. It is 
clear from Figure [TD] that the H-layer is FUV-continuum dominated, whereas the molecular 
disk (and therefore the vast majority of disk mass) is Lya-dominated. Ultimately, the 
enhancement of the Lya photon density (above the intrinsic stellar value) may exceed an 
order of magnitude. The ratios eventually asymptote to a constant value deep into the disk 
where both fields are approaching the diffusive limit and therefore behaving similarly. 

5. Discussion 

Regardless of the value of e, essentially every stellar Lya photon is intercepted by the 
highly flared H-layer. Upon the first scattering the Lya field is completely isotropized and 
a fraction (< 50%)is transmitted through to the base H-layer. This is clearly seen in Fig. [9] 
(top left and bottom left panels). Viewed from within the molecular region the Lya would 
seem to originate in a diffuse blanket overlying the disk, analogous to the diffuse transport 
of sunlight on an overcast day. This is similar to the pen etrative advantage that i nterstellar 



Willacy &: Langer 



2000l ). however 



photons (bathing the disk isotropically) would have (e.g. 

in our models the stellar Lya field is typically orders of magnitude more intense. Once in 

the realm of the molecular disk, Lya transport will be controlled primarily by dust. 
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FUV-continuum photons pass unimpeded through the H-layer until they reach the 
dust irradiation surface. Until this point the field only experiences geometric inverse- square 
law diminution - in contrast to the Lya, which has been additionally processed by 
the H-layer. In this regime the Lya/FUV-continuum photon density ratio favors the 
unimpeded FUV-continuum photons. The FUV-continuum eventually impinges upon the 
dust irradiation surface, doing so at very small angles (< 0.05 degrees). The relative 
inefficacy of scattering by dust grains is exacerbated by accompanying absorption. It is 
straightforward to show with ID planar radiative transfer models that the energy density 
transmitted downward s from this dust irrad iation surface is typically of order < 1% of 



that incident upon it ( IChandrasekhar 



19601 ). Consequently, the Lya/FUV-continuum 



photon density ratio undergoes a dramatic reversal in the vicinity of the dust irradiation 
surface, resulting in a Lya-dominated field. This resurgence appears to more than 
compensate for the reduction of Lya density due to the H-layer. In the asymptotic limit 
(large optical depth) of the analytic planar model both fields are attenuated according 
to Uph oc exp [— /cr^] , where r^ is the vertical optical depth due to dust, and fc is a decay 



constant that depends only on the scattering properties iuj^g) of dust fiFlannery et al. 
19801 ). As observed, the Lya/FUV-continuum ratio approaches a constant value deep inside 
the disk (the 'attenuated region' in Figure [TOj) . At this point both Lya and FUV-continuum 
fields are largely isotropized and their transport controlled by dust alone. In effect, the 
contrasting early life-histories are imprinted into boundary conditions that determine 
the proportionality coefficients in the asymptotic regime. A simplified diagrammatical 
representation of the transport routes is shown in Figure O 

As mentioned previously, the D'Alessio disks to not include separate thermodynamic 
treatments for gas and dust. In effect they assume that the two components are thermally 
coupled. Models show that this couphn g breaks down at low densities, causing the gas 



temperature to rise above that of dust (JGlassgold et al. 



20041 ). In turn this increases the 
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local scale height of the gas while decreasing its density. Although not explored in this 
paper, the thermal decoupling will most likely increase the extent of atomic H in the upper 
layers of the disk. This is due to two factors; first, the formation rate of H2 is density 
dependent and an increase in the disk scaleheight will be accompanied by a decrease in 
volumetric density. Second, H2 formation on grains becomes increasingly inefficient at 
T > 500- ftr due to the short time that H atoms spend on grain surfaces before thermally 



escaping (ICazaux fc Tielens 



200J). 



The settling of dust affects the Lya and FUV-continuum fields in two distinct ways. 
First, the removal of dust decreases the efficacy of the H2 formation process, resulting in a 
thicker H-layer (Eqn. [2l Fig. [8]). This ensures that the H-layer processes every stellar Lya 
incident upon it. Second, the irradiation surface for FUV-continuum photons recedes into 
the disk as e decreases. The disk becomes flatter and intercepts a smaller fraction of stellar 
FUV-continuum photons. A self-consistent disk model that couples disk structure with our 
new radiative transfer results is required in order to quantify these effects more precisely. 

We conclude with a brief summary of the main results: 

1. The vertical proflle of the Lya/FUV-continuum photon density ratio is strongly 
stratified, and may vary by orders of magnitude relative to the intrinsic stellar ratio. 
The stratification is strongest in the inner disk. 

2. The radiation field in the vicinity of the disk midplane - where most of the disk 
mass resides - is relatively enhanced in Lya. This enhancement is in addition to that 
already present in the stellar spectrum. 

3. If previous radiative transfer calculations most closely resemble our FUV-continuum 
results, then the omission of resonantly scattered Lya implies that the mass-averaged 
intradisk FUV field has been systematically underestimated. This directly affects 
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estimates of photoelectric heating rates, thermal balance, and the resulting hydrostatic 
structure. 

4. The settling of dust (e < 1) lowers the dust irradiation surface and thickens the 
resonant scattering H-layer by inhibiting H2 formation. Both effects contribute to 
a more stratified Lya/FUV-continuum ratio, ultimately resulting in a greater Lya 
enhancement in the molecular disk. 

5. The presence of an optically thick, high-albedo H-layer implies that the direct 
observation of Lya scattered by the disk (> 50% of incident photons) may be possible 
in sufficiently close, low extinction systems (e.g. TW Hya). These observations may 
place constraints on the shape and thickness of the H-layer. 

6. A Lya-dominated radiation field drives differential photo-chemic al processes in both 



the gas and grain surfaces (JBergin et al. 



2003 



Fogel et al 



2011 



). Numerous pairs 



of closely-related mol ecules exhibit a similar d ifferential response to the spectral 



form of the UV field fivan Dishoeck et al. 



20061 ). In view of the strongly stratified 



Lya/FUV-continuum ratio one might also expect the photochemical environment to 
be similarly stratified. 
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Monte Carlo radiative transfer 



Markov Chain Monte Carlo simulation is an intuitive and versatile technique for 
numerically solving high-dimensional, geometrically complex transport problems. Although 
the method permits many degrees of abstraction, the most intuitive implementation simply 
involves following the stochastic life-histories of many discrete particles (in our case, 
photons or photon packets) as the y travel th r ough a mediurn, exp e riencing interactions 



that are modeled probab i 



Bjorkman fc Wood 



2001 



isticallv (IWitt 



Gordon et al 



1977; 



2001 



Audic fc Frisch 



Whitney et al. 



1993 



Code fc Whitney 



1995 



20031 ) . If we consider the 



basic problem of photons being absorbed and scattered then it is primarily the inclusion of 
scattering that complicates the calculation, and thus motivates the use of the Monte Carlo 
method. Additionally, the Monte Carlo method places essentially no constraints on the 
geometry involved, provided a meaningful spatial discretization can be obtained. 



Witt k GordonI (119961 ). 



Our Monte Carlo simulation method is conceptually similar to 
The reader is referred to this paper for a description of the fundamental steps, which we 
briefly enumerate here. 



1. The i stellar photon packet (from a total of Nph) starts at r = 2; = with a 
luminosity Wi = l/Nph and direction cos a = i/Nph, where a is its elevation angle 
above the midplane. This is recast into a cartesian direction vector k = {k^, ky, k^) 
where k^. = cos a, ky = 0, and k^ = sin a. The star is an isotropic light source with 
total luminosity of L^ = ^^ Wi = 1. 

2. Randomly sampling Beer's Law generates the optical depth to the next scattering 
event, Tgcat = — ln(l — p), where p is a uniformly distributed random number in the 
range [0, 1]. 



3. The photon trajectory is propagated through the medium until it has traversed an 
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optical distance Tgcat = J /^ds, where the opacity k = cr^^"*rajj for FUV-continuum 
photons and k, = crJ^'^t^H ~'~ "^Lva^^-^) ^*^^ ^J(^ photons. The photon is transported 
from cell- face to cell-face (see Appendix [B]), each time traversing some distance As 
and optical distance kAs, until the scattering location is overstepped, at which point 
the photon is regressed back to the exact scattering location. The optical properties 
within a cell are considered to be homogeneous. Although the photon is traveling in 
3D space, at each face-intercept its coordinates and direction vector are rotated back 
into the r — z plane. Relevant information (e.g. packet luminosity Wi) is deposited 
as the photon passes through each cell. There is a continual attrition of photons due 
to absorption, causing a decrease in the weight-factor Wi — )■ exp[((X' — l)ijj~^Tscat]Wi, 
where u is the albedo. 

4. The scattered photon direction is found by performing these steps; transforming 
into a frame aligned with the photon direction vector, finding the scattering angle 
by randomly sampling the scattering phase function, and finally transforming 
back into the coordinate frame of the disk. This coordinate transformation is 
an application of Euler angles. In the case of Lya, the type of scattering must 
first be resolved. If a uniformly distributed random variate p G [0, 1] satisfies 

p < o"Ly^^(II)/(cr^^"*?2j| + (TLy^n(II)), then we execute resonant scattering, otherwise 
the scattering is by dust. 

5. Repeat from step 2 until photon leaves domain. 

The Henyey-Greenstein phase function used for dust scattering has a number of 
convenient features that make it attractive for use in Monte Carlo simulations. Not only 
can its probability distribution (Eqn [T]) be integrated over angle to produce a closed-form 
cumulative probability distribution, but the result can be algebraically inverted, permitting 
the expression of the scattering angle, fi = cos 6, in terms of a uniformly distributed random 
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number, p G [0, 1], 



2^ 



9 / -L — 9^ 



l + ^(2p 



(Al) 



We approximate Lya scattering as isotropic. This simplifies matters considerably, since 
the new scattered direction k is simply the result of randomly sampling directions over the 
unit sphere, thus avoiding the computationally expensive Euler angle calculation. 

We restrict our computational domain to 2; > i.e. the volume above the disk 
midplane. This is equivalent to asserting a mirror plane at z = 0; if any photons reach the 
disk midplane they are reflected k^ — ?■ —k^. The solution is now smooth and symmetric 
about the disk midplane, properly accounting for photons that pass vertically from one 
side of the disk to the other. When projected onto the r — z plane, straight lines typically 
appear as hyperbolae. Photons may escape by reaching the outer and upper edges of the 
domain. In fact, the Delaunay tessellation produces a convex domain perimeter, which in 
general is not restricted to a rectangular formo 



B. Intersection between photon path and cell walls 

Although the unstructured grid is specified in 2D, we still need to envision the cells as 
3D entities in order to use them in the Monte Carlo radiative transfer simulation. This is 
because a photon travels in 3D (it is not constrained to a single r — z plane). By employing 
an axisymmetric cylindrical framework we are asserting invariance upon rotation about 
the z axis. Once rotated, the triangular simplex is extruded into a torus with a triangular 



^The outer perimeter - or 'convex hull' - is the outline that would be formed by a taught 
string enclosing the entire set of nodes. 



-24- 

cross-section. Photon trajectories will be propagated through these volumes, requiring 
geometric tools that allow us to keep track of the intercepts between the photon path and 
the surfaces of these contiguous volumes. 

Figure [11] shows how the triangular torus is formed from the interstice of three 
intersecting cones. Computationally, the most efficient way to transport photons through 
these volumes is to jump from one cell face (the entry face) to another (the exit face). In 
order to determine these locations we need to calculate the intersections between a straight 
line (photon trajectory) and the three cones that form the surfaces of our toroidal cells 
(in 3D). The exit face is always the face with the intersection point that is closest to the 
current position of the photon (in the direction of photon propagation). 

The surface of a cone satisfies the equation, 

z ■ (x — v) = |x — v| cos6', (Bl) 

which simply states that the vector between a location in space x and the cone vertex v, 
makes an angle 6 with the cone axis (in our case the z-axis, with unit direction vector 
z). We refer to 9 as the opening angle. It is convenient to square Eqn. IBll to obtain a 
quadratic version of the cone equation. There will now be two solutions, one of which 
corresponds to the refiection cone. For example, if x lies on the cone, then so too does 
the point 2v — x that simply corresponds to the refiection of x through the vertex v. 
This formulation automatically deals with cones that are 'upside down' (i.e. have opening 
angle 9 > 90deg). Thus, the formulation is comprehensive and it only remains to select 
the solution corresponding to the appropriate cone orientation. We proceed by writing the 
squared Eqn IBll in matrix form, 

(x - v)^Af (x - v) = 0. (B2) 
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The matrix M = (zz-^ — cos 6'^ I) where / is the unit matrix. The straight hne trajectory of 
a photon can be expressed parametrically as, 

x(t) = p + tk, (B3) 

where k = {k^, ky, k^) is the unit direction vector, p is a point of origin (the current photon 
location), and t the distance along the trajectory. The intersection between the photon 
trajectory and cell wall is found by inserting Eqn IB3I into Eqn IB2I and solving for the 
distance t. Upon doing this we obtain the quadratic equation, 

C2t^ + cit + Co = 0, (B4) 

where C2 = k"^Mk, ci = k"^M (p — v), and cq = (p — v)^M(p — v). The most case is that 
C2 7^ 0, making Eqn IB4I a quadratic equation with solutions t = (— ci ± 2y^cl — C2Co)/2c2. 
If ci — C2C0 < no real solutions are possible and no intersection exists. If cf — C2C0 = 
a single (repeated root) intercept is found, which is seen to be tangent to the cone at the 
point t = — C1/C2. If cf — C2C0 > two intercepts exist at — ci ± a/c^ — C2C0/C2. There are 
additional possible situations, for example, C2 = renders Eqn IB4I linear, implying that only 
one intersection exists. In this rather improbable case the photon trajectory is parallel (but 
not coincident) with a straight line on the cone. Finally, if C2 = Ci = Cq = the photon path 
lies on the cone surface, which is exceedingly unlikely. A computational implementation 
should takes these cases into account. 

There are extremal cones that may arise from the tessellation and need to be taken 
into consideration. If the cell edge in the r — z plane is perpendicular to the z-axis, 
z ■ (A — B) = 0, then upon rotation we have a horizontal plane rather than a cone. This 
can be seen by putting 6 = 90 into Eqn IB1[ Finding the intersection of a line and a plane 
is straightforward. On the other hand, if the cell edge in the r — z plane is parallel to the 
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z-axis, for example z ■ (A — B) = |A — B|, upon rotation we have a cylinder instead of a 
cone. 
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center (jHerczeg et al. 
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rum of TW Hya expressed in Doppler units, x, relative to line 



2002) ■ Bottom - The resonant scattering Voigt profile {solid line) 
presented by a gas of H atoms with kinetic temperature T = lOOOK. The Lorentzian wings 
dominate the profile at displacements of |x| > 3. Dust opacity {dashed line) is essentially 



constant over this small wavelength range, scaling in proportion to the dust abundance e. 
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Fig. 2. — The total hydrogenic density distribution, njj, in the e = 0.01 modified D'Alessio 
disk {filled contours). A hypothetical atomic hydrogen layer ('H layer') is shown by the 
hatched region. Stellar Lya photons strike the upper surface of this layer. A fraction of 
this fiux is transmitted, emerging isotropically (shown by arrows). In contrast, stellar FUV- 
continuum photons are shown to strike the irradiation surface determined by dust opacity. 
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Fig. 3. — Optical properties of the grain ensemble composed of silicon and graphite compo- 
nents (ref to D'Alessio). Top - absorption cross-section. This paper is primarily concerned 
with dust opacity at the Lya wavelength (1215. 67A, denoted by vertical dotted line). Middle 
- scattering cross-section. Bottom - single-scattering albedo, u, and asymmetry parameter, 
g. At 1215.67Adust has an asymmetry parameter of (? ~ 0.6 and albedo u ~ 0.4. 
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Fig. 4. — The Henyey-Greenstein phase function evaluated for g = 0, 0.4 and 0.6 (see Eqn[T]). 
Isotropic scattering has g = and is a convenient approximation for modehng the resonant 
scattering of Lya by H atoms. The dust ensemble used in this paper has an asymmetry 
parameter of (7 ~ 0.6 and is therefore considered to be significantly forward-throwing. 
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Fig. 5. — Adaptation of the unstructured grid during the iterative H/H2 calculation in the 
e = 0.01 modified D'Alessio disk. Left - starting grid comprised of nested regular background 
arrays and the original D'Alessio 1+lD data. Total number of nodes is ~ 4 x 10'^. Middle 
- after one iteration new nodes have been inserted into regions where a change in the H2 
molecular fraction has been detected. Middle - the node distribution after two iterations. 
The region requiring refinement becomes smaller as the H/H2 transition becomes better 
defined. The H/H2 distribution typically converges to a solution after about 15 iterations, at 
which point the total number of nodes is ~ 3 x 10^. In the inner disk the spatial resolution 
of the discretization is on the order of 
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Fig. 6. — Detail of the inner regions (r ~ lAU) of the discretized e = 0.01 modified D'Alessio 
disk prior to the H/H2 calculation. Shown are the actual connections between nodes calcu- 
lated by performing a Delaunay tessellation on the node list. The resulting triangles are the 
cells through which photon packets propagate. As in Figure the regular array of nodes 
is the background grid, employing a divide-by-2 scheme to increase the resolution nearer to 
the star. Note that although these nodes are distributed in a rectangular array, the resulting 
cells generated by the Delaunay tessellation are triangular. The original data from the 1+lD 
D'Alessio disk model are identified as the vertical columns of nodes disrupting the regularity 
of the 2D background lattice. In a sense, the Delaunay tessellation heals these discordances, 
merging disparate sets of nodes. 
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Fig. 7. — Photodissociated layer of atomic H in the three modified D'Alessio disk models. 
The filled contours represent the total density, nu. The atomic H density n(H) is shown by 
the unfilled contours, each contour being separated by a factor of 10°'^. The solid contour 
corresponds to n(H) = lO^cm^. 
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Fig. 8. — The vertical distribution of atomic and molecular hydrogen at r = lOOAU in the 
modified D'Alessio disks with varying degrees of dust settling, e. Vertical column densi- 
ties of atomic hydrogen (A^(H), dashed line) and molecular hydrogen (A^(H2), solid line), 
and (scaled) volumetric atomic hydrogen density (^(H), dot-dashed line) are shown against 
the total hydrogen column, A^jj. Vertical column densities are integrated downwards from 
infinity. 
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Fig. 9. — Maps (r — z) of the photon density and anisotropy of the Lya {left-hand panels) 
and FUV-continuum {right-hand panels) radiation fields in the e = 0.01 modified D'Alessio 
disk model. The star is located at the origin. Top - photon density (units cm~^). Black 
arrows indicate the direction of the net flux. Bottom - the anisotropy of the radiation field, 
7 =1 J/(k)kc?i7 I /47rJ, where k is the unit direction vector, /(k) is the intensity, and J is 
the mean intensity. The limit 7 = 1 implies unidirectionality and 7 = implies isotropy. 
The flared feature of low anisotropy seen in the Lya panel {lower left) is attributable to (and 
coincident with) the resonantly scattering H-layer. Note also how the net Lya flux emerging 
from the base of the H-layer does so perpendicularly. 
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Fig. 10. — The ratio of Lya and FUV-continuuin photon densities as a function of vertical 
column density in the modified D'Alessio disks. Top - at r = lAU the ratio departs from 
the intrinsic stellar ratio (~6) as we descend into the disk (increasing iVjj). The asterisks 
denote the lower surface of the H-layer (c.f. Fig. [8]). Bottom - the same at r = lOOAU. The 
shaded areas provide an estimate of the Monte Carlo noise. The curves are terminated at 
the depth where the noise begins to degrade the data significantly. Interpretative text has 
been added for the r = lOOAU, e = 0.01 case. 
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Fig. 11. — In the r-z plane the nodes A, B and C are connected by the Delaunay tessellation 
to form a triangular cell AABC. The three extended lines that form the cell edges intercept 
the z-axis at points fi,f2 and ^3. Upon rotation about the z-axis these lines become the 
surfaces of three cones with vertices at Vi, V2 and ^3. The cone with vertex vi is shown 
by the dashed lines. The interstice formed by the three intersecting cones is a toroid with 
triangular cross-section. 



